Qubit-oscillator system: An analytical treatment of the ultrastrong coupling regime 
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We examine a two-level system coupled to a quantum oscillator, typically representing experi- 
ments in cavity and circuit quantum electrodynamics. We show how such a system can be treated 
analytically in the ultrastrong coupling limit, where the ratio g/Q between coupling strength and 
oscillator frequency approaches unity and goes beyond. In this regime the Jaynes-Cummings model 
is known to fail, because counter-rotating terms have to be taken into account. By using Van Vleck 
perturbation theory to higher orders in the qubit tunneling matrix element A we are able to en- 
large the regime of applicability of existing analytical treatments, including in particular also the 
finite bias case. We present a detailed discussion on the energy spectrum of the system and on the 
dynamics of the qubit for an oscillator at low temperature. We consider the coupling strength g 
to all orders, and the validity of our approach is even enhanced in the ultrastrong coupling regime. 
Looking at the Fourier spectrum of the population difference, we find that many frequencies are 
contributing to the dynamics. They are gathered into groups whose spacing depends on the qubit- 
oscillator detuning. Furthermore, the dynamics is not governed anymore by a vacuum Rabi splitting 
which scales linearly with g, but by a non-trivial dressing of the tunneling matrix element, which 
can be used to suppress specific frequencies through a variation of the coupling. 

PACS numbers: 03.67.Lx, 42.50.Pq, 85.25.Cp 



I. INTRODUCTION 



The model of a two- level system coupled to a quantized 
oscillator experiences widespread application in many 
different fields of physics. In quantum optics it describes 
the interaction of light with matter - of an atom coupled 
to the electromagnetic mode of a cavity. Most interest- 
ing in this instance is the regime of strong coupling; i.e., 
the coupling strength g between the atom and the cavity 
mode exceeds the loss rates stemming from spurious pro- 
cesses like escape through the cavity mirrors, relaxation 
to other atomic levels or into different photon modes, or 
decay due to fluctuations in the qubit control parameter 
induced by the environment. Under this condition, the 
atom and the cavity can repeatedly exchange excitations 
before decoherence takes over. The resulting Rabi oscil- 
lations have been observed experimentally and the field 
is known today as cavity quantum electrodynamics [l|, y] ■ 
But also for artificial atoms, like superconducting qubits 
[j-Q, similar setups have been realized with the cavity 
being formed by a one-dimensional transmission line res- 
onator [a, 0] or a simple LC-circuit [^ Q . In both cases 
the Rabi splitting in the qubit-oscillator spectrum could 
be detected [3,Q i while in the experiment of Johansson et 
al. [9[ coherent vacuum Rabi oscillations were observed. 
The advantages of this field, known as circuit quantum 
electrodynamics (QED), are manifold: For instance, the 
transition dipole moment of a superconducting Cooper- 
pair box can be made up to four orders of magnitude 
larger than in real atoms. Using a coplanar waveguide 
as cavity, the volume can be confined very tightly in the 
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transverse directions only limited by the qubit size, which 
can be made much smaller than the resonator wave- 
length. Thus, we can speak of a quasi-lD cavity, which 
leads to a strongly enhanced electric field [0, 0] and the 
strong coupling limit is more easily reached. In the first 
realization of Wallraff et al. [7| a coupling strength of 
g/n ^ 10~^ was observed, while in more recent experi- 
ments co uplings up to a few percent, g/n < 0.025, were 
reported jlOMl4| , reaching the upper limit possible for 
electric dipole coupling [1^, [T^ , whereas in cavity QED 
one finds typically g/fl ~ 10~^ [l|. The artificial atom 
can be placed at a fixed location in the cavity, so that 
fluctuations in the coupling strength are avoided. Fur- 
thermore, fabrication techniques known from integrated 
circuits can be used to "wire-up" the qubit cavity sys- 
tem and connect it to other circuit elements [ig. For 
investigations on the qubit-oscillator setup, the Jaynes- 
Cummings model (JCM) [131 i^ usually invoked. It re- 
lies on a rotating- wave approximation (RWA), which is 
valid for not too strong coupling g <C Ab,f2 and weak 
detuning, Ab ~ ^, where the qubit transition frequency 
Ab = Ve^ + A^ equals the tunneling matrix element A 
for zero static bias e. However, for certain experimen- 
tal conditions, coupling strengths of more than a few 
percent or even unity were predicted reaching the ultra- 
strong coupling regime [1^, [la, [iM [3 ■ ^or those strong 
couplings, the application of a rotating-wavc approxima- 
tion and thus the JCM is not justified anymore. For 
instance, quite recently an experiment by Niemczyk et 
al. [201 could show the failure of the JCM for a Joseph- 
son flux-qubit placed inside the center conductor of an 
inhomogeneous transmission-line resonator. Also for a 
fiux-qubit coupled to an LC-circuit, the break-down of 
the rotating wave approximation has been demonstrated 
experimentally [2l| and the ultrastrong coupling regime 
seems to be in close reach ^3. While in the JCM the 



ground state of the qubit-oscillator system consists of a 
product of the qubit's ground state and the oscillator's 
vacuum state, an inclusion of the counter-rotating terms 
leads to - depending on the coupling strength - an en- 
tangled or a squeezed vacuum state containing virtual 
photons [i^, l2J|, which under abrupt switch-off of the 
coupling arc emitted as correlated photon pairs, remind- 
ing of the dynamical Casimir effect ,19, 24, 2f|. Such an 
adiabatic manipulation has been recently realized exper- 
imentally for intersubband cavity polaritons in semicon- 
ducting quantum wells |24| . In this experiment and also 
in [2g a dimensionless coupling strength of about 10% 
has been reached. Furthermore, ultrastrong coupling has 
been predicted for qubits coupled to nanomechanical res- 
onators [26l |. 

Theories examining the qubit-oscillator system going be- 
yond the RWA are at hand: The adiabatic approximation 
(see [2Q| and references therein) relies on a polaron trans- 
formation and is derived under the assumption fi ^ Ab- 
It fails to return the limit of zero coupling g -> 0, where 
the JCM works well. An improvement to this theory 
is given by the generalized rotating-wave approximation 
(GRWA) [231, which is a combination of the adiabatic 
approximation and the standard RWA and works well in 
both regimes of zero and large qubit-oscillator detuning. 
Further, it covers correctly the weak coupling limit. How- 
ever, it has not been used yet to investigate the dynamics 
of the qubit-oscillator system. The NIBA calculations by 
Nesi et al. [23| treat analytically a two-level system cou- 
pled to a harmonic oscillator to all orders in the coupling 
strength g, taking environmental influences into account. 
Zucco et al. present a theory beyond the rotating-wave 
approximation in the strong dispersive regime [29[ . From 
these works, one can learn that the simple picture of 
the qubit-oscillator energy spectrum is not given by the 
Jaynes-Cummings ladder anymore, where pairs of energy 
levels which arc degenerate for g = are split by 2(7-y/j, 
with j denoting the higher oscillator level being involved. 
However, all these theories are derived for an unbiased 
qubit (e = 0) or in the terminology of cavity and circuit 
QED, for a qubit operated at the degeneracy point or 
sweet spot. While this situation is usually encountered 
for real atoms in cavity QED, it is quite straightforward 
to vary the static bias e of superconducting qubits by 
an external control parameter like the gate voltage ap- 
plied to a Cooper-pair box or the magnetic flux acting 
on a Josephson junction. Indeed, such a detuning from 
the degeneracy point is performed in spectroscopic mea- 
surements of the qubit-oscillator system, see e.g . [3, [21| , 
or in a current-based read-out of the qubit |30[ . There- 
fore, theories are necessary which treat the biased qubit- 
oscillator system in the ultrastrong coupling limit. In 
|3ll |33 | this is done for a qubit coupled to a linear or non- 
linear oscillator, respectively, up to second order in the 
coupling strength g. Higher-order effects like the Bloch- 
Siegert shift of the qubit dynamics could be observed. 
Brito et al. used in [33 a slightly changed polaron trans- 
formation on the qubit-oscillator model and obtained by 



truncating the displaced harmonic oscillator to its first 
excited level an effective four-level model. Quite recently, 
the adiabatic approximation for a high-frequency oscilla- 
tor was reviewed for a biased system [23| . Furthermore, 
the opposite regime of a high-frequency qubit has been 
examined there. 

In this work, we present a theory which takes the static 
bias of the qubit into account and treats the qubit- 
oscillator system to all orders in the coupling strength. 
We consider the qubit tunneling matrix element A as a 
small perturbation. For zero static bias, our approach 
can be seen as an extension of the adiabatic approxima- 
tion by taking into account higher order terms of A using 
Van Vleck perturbation theory (VVP). We do not only 
examine the energy levels of the system but also calculate 
corrections to the displaced qubit-oscillator states, which 
we obtain using a polaron transformation on the unper- 
turbed (A = 0) case. Unlike in the adiabatic approxima- 
tion discussed in [23, we take the qubit's static bias into 
account while identifying degenerate subspaces, thereby 
adjusting the renormalized frequency already in the first 
order approach. Our results work very well for nega- 
tive detuning (Ab < ri) for the whole range of coupling 
strength and even exceeds in accuracy results obtained 
from the GRWA for e = 0. For not too weak coupling 
g/fl > 0.5 and/or finite static bias, it agrees with numer- 
ical results even for the resonant case Ab = i7 or positive 
detuning Ab > fi. With these observations we believe 
we can close the gaps which cannot be treated by the 
Jaynes-Cummings model or the GRWA. 
With our investigations we enter a new physical regime: 
the splitting between the energy levels does not scale lin- 
early in g anymore but depends through a dressing by 
Laguerre polynomials on the coupling strength. This de- 
pendence allows for a suppression of individual frequency 
contributions to the dynamics. We further discover that 
even at low temperatures several frequencies come into 
play, while the JC dynamics is usually governed by two 
main oscillations. 

The outline of this work is as follows: After introduc- 
ing the Hamiltonian of the qubit-oscillator system in 
Sec. Ill Al we explain how it can be approximately diag- 
onalized by a combination of displaced oscillator states 
and VVP. The resulting eigcnstates and eigcnenergies are 
given in Sec. IIIBI being valid for the zero and nonzero 
bias case. For both situations, we examine the energy 
spectrum in detail in Sec. IIIIl comparing the different 
approaches to numerical calculations. In Sec. IIVI we 
concentrate on the dynamics; i.e., we determine the time 
evolution of the population difference of the two- level sys- 
tem and test the adiabatic approximation and VVP again 
against numerics. We conclude our discussion in Sec. |V| 



II. DIAGONALIZATION OF THE 
QUBIT-OSCILLATOR HAMILTONIAN 

A. The two-level-oscillator Hamiltonian 

The predominant model to describe the interaction be- 
tween an atom and the field of a cavity is the two-level- 
oscillator Hamiltonian, see, e.g., |34l |. 



H = i^TLS + H\ 
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(1) 



The atom is described as a simple two-level system 
(TLS), 



H- 



TLS 



\{e(T^ 4- Acr^;), 



(2) 



where we use as basis the so-called localized states, which 
are eigenstates of the Cz Pauli matrix, Uz | t / i) = i I t 
/ \). Tunnelin g b etween the two states is taken into 
account by Acr^ [421 , and e describes a possible static bias 
of the TLS. In cavity QED setups one typically finds the 
situation of zero static bias, while in circuit QED e can 
be controlled in situ. The atom is connected to the field 
of the cavity via a dipole coupling, which is expressed by 



iJint ^hga^b^ + b). 



(3) 



The coupling strength is given by g, while b'' and b arc 
the raising and lowering operators of the field. As usual, 
we assume that this field can be expressed by a single 
harmonic oscillator mode of frequency fl, 



Hn, 



hnb^b, 



(4) 



where we neglected the zero-point energy. Despite 
its simplicity, this Hamiltonian cannot be diagonalized 
analytically, and several approximation schemes have 
been developed. The most famous one is the Jaynes- 
Cummings model [17|, which neglects "energy non- 
conserving" or counter-rotating terms, and is restricted 
to relatively weak coupling strengths g <t^ Ab,ri, where 
Ab = Ve^ + A^. and to systems close to resonance, 
Ab ~ 57. A natural extension to the Jaynes-Cummings 
model ( JCM) is given in [3l[ , where the counter-rotating 
terms in the Hamiltonian ([1]) are taken into account by 
using VVP to second order in the qubit-oscillator cou- 
pling. This method thus works also for intermediate cou- 
pling strengths and biased qubits and is able to explain 
effects which go beyond the capabilities of the JCM like 
the Bloch-Sicgcrt shift recently measured in 2l|. An 



approach which goes beyond the restriction of weak cou- 
pling is the "adiabatic approximation in the displaced 
oscillator basis" , see [2g and references therein. It is de- 
rived for the limit f2 ^ Ab and relies on a separation of 
timescales: In order to calculate the fast dynamics of the 
oscillator (fast compared to the qubit), the part coming 
from the TLS in Eq. ([T]) is neglected, so that one gets an 
effective Hamiltonian for the oscillator reading 



hga^{b^ +b) + hnb^b. 



(5) 



Thus, depending on the state of the qubit the oscillator 
is displaced in opposite directions, while not changing its 
energy for a fixed oscillator quantum j, as its eigenener- 
gies are given by hjVt — hg^/Q? [2^ ■ By reintroducing the 
qubit contribution this degeneracy is lifted. However, as 
long as Ab ^ O, the doublet structure is conserved. For 
an unbiased system, as done in [26| , the condition trans- 
lates to A ^ J7 and the tunneling matrix element A can 
be treated as a small perturbation, in the end leading 
to an effective Hamiltonian consisting of 2 by 2 blocks, 
with a renormalized frequency on the off-diagonal. As 
this special case is included in our calculation, we will 
describe it in more detail below. Furthermore, the con- 
trary regime of a high-frequency qubit Ab ^ 0, has been 
treated in [23 analytically for certain special cases. This 
situation is also partly contained in our formalism. 



B. Eigenenergies and eigenstates 

In the following, we demonstrate how the full Hamil- 
tonian H can be diagonalized pcrturbatively to second 
order in A. For a vanishing tunneling element, A = 0, 
the polaron-likc transformation 



U 



„g(6-&t)CT,/n 



(6) 



brings H into a diagonal form [43|. Its eigenstates are 

It / i,j) = U\\ / i,j), where | t / i,i) arc the eigen- 
states of the qubit-oscillator system for A = and g = 0. 
For detailed expressions see Eqs. (|A1[) and (|A2p . They 
correspond to the displaced oscillator states used in l2g , 
where the displacement depends on the qubit state. The 
eigenvalues are 









(7) 



For finite A, the perturbative matrix elements become 



2 ^ 



{i,j\AaM,j') 
A[sign(/-j)]l^'--'lsl^;3,^,j(a), (8) 



with 



S'(a) = a'^ 



U + i) 



^(a)e- 



(9) 



and a = {2g/il.)'^. This dressing by Laguerre polynomials 
becomes in the high photon limit, j — > (X), and for finite 
I a dressing by Bcssel functions, just like in the case of a 
classically driven TLS ^^^. 

For A = and e = /fi, the unperturbed eigen- 
states \i,j) and It, i + arc degenerate, so that we 
can identify a twofold degenerate subspace in the com- 
plete Hilbcrt space of the problem [4J] . By using VVP 



[40|, we can determine an effective Hamiltonian Hcff ~ 
exp(i5')i7 exp(— iiS") for the perturbed system consisting 
of 2 by 2 blocks of the shape 



^L- 



/i,(2) 



2^i 



i^f.J+l 



(10) 



where we calculate the transformation matrix S to second 
order in A [43 and define the diagonal corrections as 



for < 9 '■ < IT. In Appendix \^ the transformation is 
calculated to second order in A and applied to the effec- 
tive states. By this we have all information we need to 
calculate the dynamics of the qubit-oscillator system. 
Van Vleck perturbation theory yields good approximate 
results as long as the matrix elements connecting differ- 
ent non-degenerate subspaces with each other are much 
smaller then the energetical distance between those sub- 
In our case this means 
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34 
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(11) 



Notice that for zero bias, e = 0, the degenerate subspace 
consists of oscillator states with equal quantum number 
j. If one neglects the second-order corrections e*-^-* the 
effective Hamiltonian reduces to the one obtained within 
the "adiabatic approximation" in [26|, see Eq. (9) there]. 
Thus, our approach automatically also includes the adi- 
abatic approximation. In [26[ only the zero bias case is 
considered; here we extend the adiabatic approximation 
to finite bias disregarding the second order correction e'^'^^ 
in Eq. (jTU]). In [231, ^ finite bias £ is considered in the 
parameter regime where eigenstates with same oscillator 
quanta j remain quasidegenerate, so that the tunneling 
matrix element of a subspace remains dressed by a L^ 
Laguerre polynomial. This is a valid approximation in 
the case that Q, ^ Ab- On the contrary, when e > fl and 
therefore also Ab > fi, a dressing by higher-order La- 
guerre polynomials occurs even in first order in A. The 
eigenenergies of Eq. ([TU| are 



Et,, = f^ 
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with the dressed oscillation frequency 
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Notice that the quantum number j corresponds to a mix- 
ture of the oscillator levels j and I. Only for e = 
this mixing vanishes. We obtain the eigenstates of H 
by |*±.j) = cxp(-iS')|$^°'j) with the eigenstates of ^ 
given by 



l$L"',) = -sin:^i;,j)- 
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n(Af')cos:^|t,j+0, 
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sign(Af jsin^|t,i + 0, 



and the mixing angle 
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(16) 
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(17) 



We will discuss the validity of our approach for the dif- 
ferent cases below. 



III. ENERGY SPECTRUM IN THE 
ULTRASTRONG COUPLING REGIME 

In this section, we examine the energy spectrum of the 
qubit-oscillator system as obtained from Eq. ([T^ and 
compare it to results found by exact numerical diago- 
nalization. We check its robustness for variable coupling 
strength g and detuning d = Ab — fi between the qubit 
energy splitting and oscillator frequency. 



A. Zero static bias £ = 

First, we concentrate on the regime of zero static bias. 
This is the usual case in cavity QED, where the JCM 
is applied. The JCM is known to work well for weak 
qubit-oscillator coupling {g/fl ^ 1) and small detuning 
between the two devices. As already predicted in j3l| . 
higher-order corrections have to be taken into account for 
stronger coupling. For the case of ultrastrong coupling, 
we will find that the situation changes dramatically. The 
energies predicted by the JCM read 



E. 



JCM 

2j + l,2j+2 



.4)^Tiy(A- 



n)^ + 4(j + l).g2 



(18) 

with the ground state energy Eq'~^^ = —hA/2. Equation 
([T2I) for the Van Vleck eigenenergies perturbativc in A, 
simplifies further for e = 0: 



Etj = ^ 
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(19) 

The semi-infinite sum in the above expression converges, 
and we show in Appendix[B]analytical expressions for the 
first four energy levels. Furthermore, we can compare our 
results to the generalized rotating-wave approximation 
(GRWA) [23|. In this approach, the total Hamiltonian 
Eq. (HI) is expressed in the displaced basis states of the 
adiabatic approximation. It is then in this representa- 
tion, that the rotating-wave approximation is performed 
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FIG. 1: Energy levels against detuning S = A — $1 for e/il = 0, 
g/Q = 0.1. Our VVP solution is compared to the GRWA and 
the JCM. The latter two agree well with numerical calcula- 
tions for the whole detuning range (not shown), while VVP 
yields only reliable results for negative detuning, A < fi. 



and counter-rotating terms are neglected. Thus, the 
GRWA uses the advantages of the adiabatic approxima- 
tion, namely its ability to go to strong coupling strengths 
and to treat detuned systems, and also gives reliable re- 
sults in the weak coupling regime of the JCM. A deriva- 
tion of the GRWA eigenenergies can be found in Ap- 
pendix [cl 



1. Energy levels against detuning 

In Figs. [T]-|3]we examine the energy levels against the 
qubit-oscillator detuning (5 = A — 51 at fixed couplings, 
g/fl — 0.1, 0.5, 1.0 and 1.5, respectively. 
For a weak coupling of g/il = 0.1, we compare VVP to 
the GRWA and the JCM. Both are known to work well 
in this regime. We find that VVP gives only valid results 
for negative detuning, A < 51. This was expected as it 
relies on a pcrturbative approach in A, and we know al- 
ready from the adiabatic approximation that it fails for 
A > il and simultaneously small g/fl. In this regime of 
weak coupling, the JCM or GRWA are clearly preferable 
to our method. 

For an intermediate coupling strength, the same discus- 
sion is presented in Fig. [2] We do not show the Jaynes- 
Cummings energy levels in this regime anymore, because 
they fail completely to return the correct energy spec- 
trum. Instead we compare to a numerical diagonaliza- 
tion of the Hamiltonian. Van Vleck perturbation theory 
and the GRWA yield good results for negative detuning 
5 < 0, but also at resonance, A = 51, they agree rela- 
tively well with the numerics. At positive detuning both 
deviate strongly from the exact solution. 
With a coupling strength of g/il = 1.0 in Fig. [31 we are 
already deep in the ultrastrong coupling regime. Those 
high values have not been observed experimentally yet. 



FIG. 2: Energy levels against detuning 5 = A — $1 for e/Q — 0, 
g/^ — 0.5. The JCM fails already completely for such a 
coupling strength (not shown). We compare VVP and the 
GRWA against numerical calculations. Both agree well with 
the numerics for negative detuning and even at resonance. 
For stronger positive detuning they both fail and strongest 
deviations can be seen for the lower energy levels. 
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FIG. 3: Energy levels against detuning 5 = A — 51 for e/Q — 0, 
g/fl = 1.0. We compare VVP, the adiabatic approximation 
and GRWA against a numerical calculation. For a negative 
detuning all three approaches agree very well with the exact 
numerics. However, for zero and positive detuning deviations 
occur. In particular, the ground level and the first excited 
level are not described correctly by the adiabatic approxima- 
tion and the GRWA for strong positive detuning, while VVP 
yields good results. 



They are, however, predicted to be realizable |15| . For 
negative detuning, GRWA and VVP show a good agree- 
ment with the numerics. However, approaching zero de- 
tuning or going beyond to positive one, the GRWA fails 
in particular for the two lowest states, which will turn 
out to be important for the calculations of the dynamics. 
In order to explain this failure, we also show in Fig. [3] 
the adiabatic approximation. As pointed out, the GRWA 
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FIG. 4: Energy levels against detuning. Same as in Fig. (3] 
but for a coupling strength of g/Q — 1.5. Adiabatic approxi- 
mation and GRWA fail for positive detuning, while VVP gives 
the first four energy levels correctly even up to a detuning of 
S/Q = 2.0. And also for the higher energy levels it yields 
good results beyond the resonant case. 



is a combination of the ordinary RWA, and thus works 
well for weak coupling, and of the adiabatic approxima- 
tion, which works very well for strong negative detuning, 
ri ^ A, for all values of the coupling. At resonance or 
at positive detuning, the adiabatic approximation shows 
deviations from the exact solution for a coupling strength 
g/il = 1.0. This coupling strength is, however, already 
too strong to be treated correctly by the RWA. Thus, 
we arc in a kind of intermediate regime, which is also 
not covered by the GRWA, but can be important in ex- 
perimental applications. On the contrary, VVP shows an 
exact agreement with the numerical data for negative de- 
tuning and even up to exact resonance. Only for positive 
detuning, deviations start to occur. 

This becomes even more prominent for stronger coupling 
strengths, like g/fl = 1.5 in Fig. JH While the adiabatic 
approximation and also the GRWA fail for positive de- 
tuning, VVP agrees surprisingly well with the numerical 
results up to ^ = 2.0 for the first four energy levels; i.e., 
we have A/ft = 3.0. Also for the higher levels we still 
find a good agreement for not too strong positive detun- 
ing. This improvement is due to the fact that VVP also 
takes into account connections between non-degenerate 
subspaccs and therefore higher-order corrections in the 
dressed tunneling matrix element. 



2. Energy levels against coupling strength 

In Figs. |S]-|7]we investigate now the first eight energy 
levels against the coupling strength g/ft for three differ- 
ent values of the detuning. 

All three approaches, the adiabatic approximation, the 
GRWA and VVP, show very good agreement with the 
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FIG. 5: Energy levels against coupling strength g/fl for nega- 
tive detuning (5/f2 = —0.5). Numerical results are compared 
with the adiabatic approximation, GRWA and VVP. All three 
approaches show only slight deviations. 
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FIG. 6: Energy levels against coupling strength at resonance 
{S/n = 0). For small coupling strength, the adiabatic ap- 
proximation and VVP show small deviations from the correct 
values (see especially the higher energy levels). The GRWA 
works well in this regime. For stronger coupling strength, all 
three approaches agree well with the numerical results. 



numerical results for the whole range of g/Q for negative 
detuning S/il = —0.5 shown in Fig. \5\ 
At resonance, A/fl = 1.0, in Fig. [HI we have to distin- 
guish between different parameter regimes: For smaller 
values of the coupling, g/fl < 0.5, the adiabatic approxi- 
mation and VVP show deviations from the numerical re- 
sults apart from the ground level, as they do not take into 
account correctly the zero coupling resonance [27[, while 
the GRWA on the other hand works well. For higher cou- 
pling strengths on the other hand, VVP exhibits a slight 
improvement to the GRWA and the adiabatic approxi- 
mation for the first two energy levels, as could already 
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FIG. 7: Energy levels against coupling strength for positive 
detuning [5/Q. = 0.5). For coupling strengths with g/Q. > 
0.75, VVP exhibits the best agreement with numerical results, 
while for smaller coupling and higher energy levels, the GRWA 
should be used. 




FIG. 8: Sketch of the validity regime of VVP and GRWA for 
£ = 0. The GRWA is perferable to VVP at weak coupling, in 
particular close to resonance and positive detuning. On the 
contrary, VVP works better at strong coupling strengths. 



be seen from Figs. [3] and 21 

This improvement becomes more evident for stronger 
positive detuning, S/i^ = 0.5, as shown in Fig. [T] 
Considering the lowest two energy levels, VVP agrees 
well with the numerical results for g/H. > 0.75, while 
the adiabatic approximation and GRWA strongly devi- 
ate from the numerical results. For higher levels also 
the latter two are closer to the numerics. However, for 
weaker couplings the results from all three approaches 
are not very satisfying even for the lower energy levels, 
and the adiabatic approximation and VVP predict un- 
physical crossings, while the GRWA at least yields the 
correct weak coupling limit. 

Plotted against the coupling strength the energy levels 
exhibit some peculiarities. Most interesting is the find- 
ing that for strong coupling two adjacent energy levels 
become degenerate, so that coherent oscillations between 
them become completely suppressed. We can understand 
that by considering expression (|19p . where we find that 
two energy levels with the same index j differ only in the 
sign of the dressed oscillation frequency, which vanishes 
for large g. For the higher energy levels, degeneracies also 
occur for lower g/Q values, happening at the zeros of the 
Laguerre polynomials. These phenomena are discussed 
in more detail in [2j, [2y, [23 , and we come back to them 
when presenting the dynamics. 



3. Validity regimes 

To summarize this section we give a comparison be- 
tween VVP and the GRWA. We do not discuss the adi- 
abatic approximation and the JCM as they are included 
in VVP and the GRWA, respectively. Further, we want 
to emphasize that Fig. [5] only represents a qualitative 



sketch; the detailed behavior is more complicated: the 
validity regime of the different approaches is crucially 
dependent on the error one allows compared to numeri- 
cal solutions. Furthermore, the number of energy levels 
taken into account plays a role. For instance, in Fig. [7] 
VVP agrees very well with the numerics for the lowest 
two energy levels and g/fl w 0.75, but shows already 
stronger deviations for the 5th and 6th level. In Fig. [5] 
we took the first eight levels into account. In order to 
understand the validity regime of VVP we consider Eq. 
(fT7)l for e — 0. In this special case it becomes 



2 ^ 



< \kn\ V k^o. 



(20) 



From the definition of the dressed tunneling matrix ele- 
ment Aj , Eq. ([HI), we see that for small A/fi - i.e., 
for negative detuning - this condition is fullfilled even 
for weak coupling. However, for A > fi and weak cou- 
pling, the above condition does not hold anymore. On 
the other hand, by increasing the coupling strength VVP 
becomes even valid at strong positive detuning since the 
dressed tunneling matrix elements are exponentially sup- 
pressed. The GRWA is valid for positive detuning also 
in the case of weak coupling. For intermediate coupling 
0-5 ^ g/^ ^ 1.0 it fails for zero or positive detuning, 
while increasing the coupling strength further yields an 
improvement in this regime. This last tendency has the 
same origin as in case of VVP, namely that the neglected 
tunneling matrix elements get suppressed. As, however, 
the GRWA considers these matrix elements only to first 
order, the improvement is not as good as for VVP. 




FIG. 9: Energy levels against static bias e for g/Q, = 1.0 
at resonance A/f2 = 1.0. Van Vleck perturbation theory is 
compared to a numerical diagonalization of the Hamiltonian. 
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FIG. 10: Energy levels against coupling g/Q, for e/fi = 1.0 
and A/fi = 0.5. The adiabatic approximation and VVP agree 
almost perfectly with numerical results. Slight deviations can 
be seen for the adiabatic approximation at g/Q. — >■ 0. 



B. Finite static bias e 7^ 

In this section wc discuss the energy spectrum for the 
case of finite static bias. We compare our VVP calcula- 
tion to exact numerical diagonalization. We further show 
in certain cases calculations disregarding connections be- 
tween the different manifolds, that is second-order correc- 
tions in A, which is the natural extension of the adiabatic 
approximation to finite bias. We do not compare to the 
GRWA, as it exists so far just for the zero bias case. To 
start, we show in Fig. IH] the energy levels against the 
static bias for a coupling strength of g/Vt = 1.0 and no 
detuning in the zero bias case (A = 57). For such a cou- 
pling strength, we find a very good agreement between 
our VVP calculations and numerically obtained results. 
Most remarkably, this agreement holds even away from 
the resonant points, e — ID,, for which our approximation 
has been performed. We also checked the effect on the 
spectrum when neglecting the second-order corrections 
in A. The qualitative behavior remains the same; how- 
ever, quantitative deviations occur (not shown in Fig. 
E]). For negative detuning, A < fl, the agreement be- 
tween analytical and numerical results is even enhanced, 
while for positive detuning up to A/fl = 1.5 only slight 
deviations occur. The accuracy of VVP diminishes en- 
tering the weak coupling regime, as we could already ob- 
serve for the zero static bias case and we will show in the 
following. Before, we want to consider some general fea- 
tures of the spectrum at nonzero static bias. We already 
pointed out while identifying the degenerate subspaces in 
Eq. ((To)) that for e ^ Ifl with I ^ certain unperturbed 
energy levels have no degenerate partner. Without loss 
of generality, we assume Z > 0, that means that the first 
I energy levels corresponding to a spin-up state have no 
degenerate partner and their energy is simply given by 
E°.^ - f £^^] with j = 0, 1, 2, ...I - 1. Of course, also the 
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FIG. 11: Energy levels against coupling g/n for e/fi = 1.0 
and A/f2 = 1.0. Van Vleck perturbation theory is still valid 
compared to numerical results, while the adiabatic approxi- 
mation fails specifically for weak coupling strengths. 



corresponding effective eigenstates are simply Itij): smd 
we cannot observe avoided crossings or a superposition 
of states. For instance, in Fig. |9] at e/ft = 1, we ob- 
serve the lowest energy level being without partner, while 
the higher ones form avoided crossings with the adjacent 
level. For e/fl = 2, the two lowest levels are "free", etc. 
In Figures [TOl [TT] and [121 we present the dependence 
of the energy spectrum on the coupling strength g/fl 
for the case of e/fl = 1.0 and A/fl = 0.5, A/0 = 1.0 
and A/n = 1.5, respectively. Just like in the zero static 
bias case, VVP yields best results for A/Q < 1, because 
there the condition for a perturbative approach is most 
satisfied. Also, the extended adiabatic approach yields 
very convincing results, only for g/D, -^ one can no- 
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FIG. 12: Energy levels against coupling g/Q, for e/fi = 1.0 
and A/r2 = 1.5. In this regime, also VVP shows deviations 
from the numerical results for g/Q, < 0.75 especially for the 
higher energy levels. It agrees well for stronger coupling. 




FIG. 13: Sketch of the validity regime of VVP and of the 
adiabatic approximation for e — 1.0. For positive detuning 
and simultaneously weak coupling both approaches fail. For 
stronger coupling VVP yields an improvement to the adia- 
batic approximation. 



tice slight deviations. For A/il = 1.0 in Fig. [TTl VVP 
still shows almost exact agreement with the numerical re- 
sults, whereas the adiabatic approximation fails for weak 
coupling. This failure of the latter becomes more evident 
going to positive detuning like A/fl = 1.5 in Fig. [T^l 
But there also the VVP exhibits strong deviations for 
coupling strengths g/il < 0.75. 

Figure [T^ summarizes these observations in a qualitative 
sketch of the validity regimes. Thereby VVP excels the 
adiabatic approximation as it considers also second-order 
corrections in the matrix elements connecting different 
doublets. 
We also tested for static bias values being no multiples 



FIG. 14: Energy levels against coupling g/Q for e/fi = 3.0 and 
A/Q — 1.5. The three lowest energy levels have no degenerate 
partner. Despite the high value of A, VVP still gives reliable 
results, while the adiabatic approximation differs from the 
numerical values even for the low energy levels. 



of fl and found a confirmation of the above findings. For 
stronger static bias, VVP describes the lower energy lev- 
els even better for positive detuning, see, e.g., the case 
e/fl = 3.0 in Fig. [M) Here, the three lowest energy lev- 
els are without degenerate partner and therefore can be 
described by the corrected unperturbed energy. The in- 
fluence of the mixing to other energy levels is less strong. 



IV. DYNAMICS OF THE QUBIT IN THE 
ULTRASTRONG COUPLING REGIME 

We are interested in determining the population differ- 
ence between the two qubit states; i.e., we calculate 

{a,{t)) = TTTLsW.Predit)} = 2(t |pred(t)| t) - 1, (21) 

where p^cdit) is obtained after tracing out the oscillator 
degrees of freedom from the qubit-oscillator density op- 
erator p. The matrix elements of the latter read in the 
system's energy eigcnbasis {|$a={±j})} 



Pajit) = ($„|p(i)|$^) = p„^(0)e-'"°-* 



(22) 



As starting conditions, we assume the qubit and the os- 
cillator to be uncoupled for t < 0, and the first to be 
prepared in the spin-up state, with the oscillator being 
in thermal equilibrium: 



p(o) = it>(ti®E|^ 



-hfijCi 



J)(J|, 



(23) 



where Z is the partition function of the harmonic oscilla- 
tor and /3 the inverse temperature. In the following, we 
will assume H^n = 10, which corresponds for oscillator 
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FIG. 15: Population difference for zero static bias. Further 
parameters are A/f2 = 0.5, hpQ = 10 and g/Sl = 1.0. The 
adiabatic approximation and VVP are compared to numerical 
results. The first one only covers the longscale dynamics, 
while VVP also returns the fast oscillations. With increasing 
time small differences between numerical results and VVP 
become more pronounced. 



FH 



dt{(Tz{t)) cos(i/t), 



(24) 



respectively. Concerning the population difference, we 
see a relatively good agreement between the numerical 
calculation and VVP for short timescales. In particu- 
lar, VVP also correctly returns the small overlaid oscil- 
lations. For longer timescales, the two curves get out 
of phase. The adiabatic approximation only can re- 
produce the coarse-grained dynamics. The fast oscilla- 
tions are completely missed. To understand this better, 
we turn our attention to the Fourier transform in Fig. 
1161 There, we find several groups of frequencies located 
around i^/il = 0, i^/il = 1.0, i^/il = 2.0 and i^/n = 3.0. 
This can be explained by considering the transition fre- 
quencies in more detail. We have from Eq. p3 
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frequencies in the GHz regime to experiments performed 
at several mK. At those low temperatures, mainly the 
lower oscillator energy levels are of importance. The dy- 
namics for higher oscillator occupation numbers at zero 
static bias has been investigated in |26| . 
The transition frequencies are defined as uja-y = (Ea — 
Ej)/h, where Ea stands either for Ez^ij in case of two- 
fold degenerate subspaces or E9,. =p l^tA i ^'^'^ '^^^^~ 
dimensional subspaces. We further can distinguish be- 
tween two different timescales: large oscillatory contri- 
butions are resulting from different oscillator quanta j, 
while the difference in dressed oscillation frequencies i7' 
acts on a much longer timescale and its contribution van- 
ishes for large coupling strengths 5/ri. 
In the following subsections we will investigate the dy- 
namics for the unbiased and biased case. Again, we will 
compare exact numerical results to VVP and the adia- 
batic approximation. Apart from the energy levels, also 
the eigenstates become now of importance. In particu- 
lar, we will find that away from the condition e = Ifl, the 
higher-order corrections are crucial to give the correct 
dynamics. 



A. Dynamics for zero static bias e = 

For zero static bias, we first examine a regime where we 
expect our approximation to work well. We thus consider 
a not too strong tunneling matrix element, A/il = 0.5 
and a coupling strength oi g/ft ~ 1.0. Figures [T51 and 
[TSlshow the population difference {<Jz{t)) and its Fourier 
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second-order corrections. For zero bias, e = 0, the in- 
dex / vanishes. The term {k — j)Vt determines to which 
group of peaks a frequency belongs and 17" its relative 
position within this group. The latter has A as an up- 
per bound, so that the range over which the peaks are 
spread within a group increases with A. The dynam- 
ics is dominated by the peaks belonging to transitions 
between the same subspace k — j — 0, while the next 
group with k — j ~ 1 yields already faster oscillations. 
To each group belong theoretically infinite many peaks. 
However, under the low temperature assumption only 
those with a small oscillator number play a role. For 
the used parameter regime, the adiabatic approximation 
does not take into account the connections between dif- 
ferent manifolds. It therefore covers only the first group 
of peaks with k — j = 0, providing the long-scale dynam- 
ics. For e = 0, the dominating frequencies in this first 
group are given by n° = |Ae-"/2|, n° = |A(1 -a)e-"/2| 
and n'^ = |AL^(a)e-"/2|, where fJg and fi^ coincide. A 
small peak at fl^ = |AL3(a)e~"/2| can also be seen. No- 
tice that for certain coupling strengths some peaks van- 
ish; like, for example, choosing a coupling strength of 
g/fl = 0.5 makes the peak at fll vanish completely, in- 
dependently of A, and the ilp and ilj peaks split. The 
JCM yields two oscillation peaks determined by the Rabi 
splitting and fails completely to give the correct dynam- 
ics, see the left-hand graph in Fig. [T51 
Now, we proceed to an even stronger coupling, g/Q = 2.0, 
where we also expect the adiabatic approximation to 
work better. From Fig. [5] we noticed that at such a 
coupling strength the lowest energy levels are degener- 
ate within a subspace. Only for oscillator numbers like 
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FIG. 16: Fourier transform of the population difference in Fig. 1151 The left-hand graph shows the whole frequency range. 
The lowest frequency peaks originate from transitions between levels of a degenerate subspace and are determined through the 
dressed oscillation frequency $1°. Numerical calculations and VVP predict group of peaks located around u/Q — 0, 1.0, 2.0, 3.0. 
The first group at v/fl = is shown in the middle graph. One can identify frequencies Qq and JI2, which fall together, and f2?. 
The small peak comes from the frequency ^3. This first group of peaks is also covered by the adiabatic approximation. The 
other groups come from transitions between different manifolds. The adiabatic approximation does not take them into account, 
while VVP does. A blow-up of the peaks coming from transitions between neighboring manifolds is given in the right-hand 
graph. In the left-hand graph additionally the Jaynes-Cummings peaks are shown, which, however, fail completely. 



. — . — . - Adiabatic Approx.- 

VVP 

Mumerical 



V 




is a clear indication that even at low temperatures higher 
oscillator quanta are involved due to the large coupling 
strength. Also frequencies coming from transitions be- 
tween the energy levels from neighboring manifolds are 
shown enlarged in Fig. 1181 The adiabatic approximation 
and VVP can cover the main structure of the peaks in- 
volved there, vifhile the former shows stronger deviations. 
If we go to higher values A/H. > 1, the peaks in the 
individual groups become more spread out in frequency 
space, and for the population difference dephasing al- 
ready occurs at a shorter timescale. For A/i7 = 1, at 
least VVP yields still acceptable results in Fourier space 
but gets fast out of phase for the population difference. 



FIG. 17: Population difference for zero static bias. Same 
parameters as in Fig. [15] but for a coupling strength of g/il = 
2.0. Both the adiabatic approximation and VVP agree well 
with the numerics, but show slight dephasing on a longer 
timescale. 



j = 3, we see that a small splitting arises. This splitting 
becomes larger for higher levels. Thus, only this and 
higher manifolds can give significant contributions to the 
long time dynamics; that is, they can yield low frequency 
peaks. Also the adiabatic approximation is expected to 
work better for such strong couplings [26[. And indeed 
by looking at Figs. [T7] and [TH] we notice that both the 
adiabatic approximation and VVP agree quite well with 
the numerics. Especially the first group of Fourier peaks 
in Fig. [TH] is also covered almost correctly by the adia- 
batic approximation. The first manifolds we can identify 
with those peaks are the ones with j = 3 and j = 4. This 



B. Dynamics for finite static bias e 7^ 

As a first case, we consider in Fig. [19] a weakly biased 
qubit {e/il = a/O.S) being at resonance with the oscilla- 
tor (Ab = rj). For a coupling strength oi g/fl = 1.0, we 
find a good agreement between the numerics and VVP. 
The adiabatic approximation, however, conveys a slightly 
different picture: Looking at the time evolution it re- 
veals collapse and rebirth of oscillations after a certain 
interval. This feature does not survive for the exact dy- 
namics. Like in the unbiased case, the adiabatic approx- 
imation gives only the first group of frequencies between 
the quasidegenerate subspaces, and thus yields a wrong 
picture of the dynamics. In order to cover the higher 
frequency groups, we need again to go to higher-order 
corrections by using VVP. 

For the derivation of our results we assumed that e is a 
multiple of the oscillator frequency fi, e = /fi. In this 
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FIG. 18: Fourier spectrum of the population difference in Fig. 1171 In the left-hand graph a large frequency range is covered. 
Peaks are located around v/n = 0, 1.0, 2.0, 3.0 etc. Even the adiabatic approximation exhibits the higher frequencies. The 
upper right-hand graph shows the first group close to ly/Q — 0. The two main peaks come from Q3 and Q4 and higher degenerate 
manifolds. Frequencies from lower manifolds contribute to the peak at zero. The adiabatic approximation and VVP agree well 
with the numerics. The lower right-hand graph shows the second group of peaks around u/Q, — 1.0. This group is also predicted 
by the adiabatic approximation and VVP, but they do not fully return the detailed structure of the numerics. Interestingly, 
there is no peak exactly at v/Cl = 1.0 indicating no nearest-neighbor transition between the low degenerate levels. 




FIG. 19; Population difference and Fourier spectrum for a biased qubit (e/Sl = -\/0.5) at resonance with the oscillator (Ab = fi) 
in the ultrastrong coupling regime ((?/n = 1.0). Concerning the time evolution VVP agrees well with numerical results. Only 
for long time weak dephasing occurs. The inset in the left-hand figure shows the adiabatic approximation only. It exhibits 
death and revival of oscillations which are not confirmed by the numerics. For the Fourier spectrum, VVP covers the various 
frequency peaks, which are gathered into groups like for the unbiased case. The adiabatic approximation only returns the first 
group. 



case we found that the levels E? 



and £'? o+; 



form a de- 
generate doublet, which dominates the long-scale dynam- 
ics through the dressed oscillations frequency J7'. For I 
being not an integer those doublets cannot be identified 
unambiguously anymore. For instance, we examine the 
case e/ft ~ 1.5 in Fig. [20l Here, it is not clear which lev- 
els should be gathered into one subspace: j and j -I- 1 or 
j and j + 2. Both the dressed oscillation frequencies fij 
and fl^ influence the longtime dynamics. In Fig. [20l we 



chose I = 2 for our approximate method. Surprisingly, 
VVP gives a very accurate picture for both the dynam- 
ics and the Fourier spectrum. For I = 1 we obtained 
the same result (not shown here). Thus, our approach 
can also treat the case of e being not a multiple of ft, 
and independent of the choice of I, VVP covers all rele- 
vant frequencies because of taking into account connec- 
tions between different manifolds. We always find pairs 
of frequencies resulting from 51^ and 51^. Those pairs 
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FIG. 20: Population difference and Fourier spectrum for e/Sl — 1.5, A/Sl = 0.5 and g/Q. = 1.0. Van Vleck perturbation theory 
is confirmed by numerical calculations, while results obtained from the adiabatic approximation deviate strongly. In Fourier 
space, we find pairs of frequency peaks coming from the two dressed oscillation frequencies ilj and Q,^. The spacings in between 
those pairs is about 0.5$!. The adiabatic approximation only returns one of those dressed frequencies in the first pair. 



are separated approximately by O.Sfi, which is the small- 
est distance between the unperturbed energy levels (only 
the single levels are separated by a larger distance). For 
a bias of e/fl = 2.5, for example, one would detect the 
same separation between the different groups of peaks. 
The adiabatic approximation extended to nonzero static 
bias fails in such a situation, as it will always only con- 
sider one of the two frequencies, which can be also seen 
by looking at the dynamics in Fig. [20l Furthermore, 
as we saw already in the unbiased case, it neglects the 
higher frequencies for intermediate coupling. 



V. CONCLUSION 

Up to now there exists no clear definition of the ultra- 
strong coupling regime. In many cases it is used to denote 
coupling strengths g/fl for which the Jaynes-Cummings 
model is not valid anymore. Consequences of this fail- 
ure can already be visualized for intermediate regimes 
like g/n « 0.1 0, iH, H ill. At such a couphng 
strength it is often sufficient to take into account the 
counter-rotating terms in the qubit-oscillator Hamilto- 
nian by treating it perturbatively to second order in g 
[3l| . For stronger coupling, like g/il, approaching unity or 
going beyond, higher orders are needed. In this work we 
presented an approach which treats the qubit-oscillator 
system to all orders in g. The price we had to pay was to 
make some restriction on the tunneling matrix element 
A and thus the qubit transition frequency Ab compared 
to the oscillator frequency fl. In detail, we followed a 
perturbative approach with respect to the dressed tun- 
neling element A-. However, since especially for strong 
coupling this dressed element becomes suppressed by a 
Gaussian, and by using VVP to include also higher or- 
ders, we could go beyond the limit A ^ fi of an adia- 



batically fast oscillator [23|, |26[. For zero bias, we com- 
pared the energy spectrum obtained by our method and 
the adiabatic approximation to the generalized rotating- 
wave approximation in [27| . For A/f2 < 1 all approaches 
agree well with numerical results for the whole coupling 
strength, while at resonance and slight positive detuning 
the GRWA was found to be preferable at weak coupling 
5 — > 0, since it returns correctly the Jaynes-Cummings 
limit. For strong coupling and small positive detuning 
VVP even showed slightly better results than the GRWA. 
We investigated in detail the dynamics of the qubit in the 
zero bias case and the ultrastrong coupling regime at low 
temperature. While the adiabatic approximation gives a 
coarse-grained picture of the time evolution of the popu- 
lation difference, VVP also covers the higher frequencies 
agreeing well with numerical results. For not too weak 
coupling our approach even gives reliable results at the 
resonance point A = J7 and slightly beyond. 
The dynamics obtained in the ultrastrong coupling 
regime is much richer than the ones predicted by the 
JCM. Instead of two dominating frequencies we found 
groups of peaks whose splitting is not linear in g any- 
more as found for the common vacuum Rabi oscillations 
but rather depends in a non-trivial way on a dressing by 
Laguerre polynomials. With the help of our analytical 
formulas we could understand this structure. The sit- 
uation reminds of the case of a classically driven TLS, 
where the resulting Rabi frequency, which for weak driv- 
ing is linear in the driving amplitude, shows a Bessel 
function like dependence in the case of extreme driving 
[33 1 . The dressing of the qubit-oscillator system by La- 
guerre polynomials allows a suppression of specific fre- 
quencies through a variation of the coupling strength g. 
Finally, we could see from the expressions ([T^ for the 
dressed oscillation frequency and (jA6p for the second- 
order eigenstates that one cannot speak of single qubit 
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or oscillator contributions anymore but has to consider a 
highly entangled system even for the ground state. 
Furthermore, we examined the situation of a biased 
qubit, which so far has not been treated analytically for 
the regime of comparable qubit and oscillator frequency 
(Ab ^ r2). An extension of the adiabatic approximation 
to the biased case was almost automatically included in 
our treatment. We showed that for situations where the 
bias is not a multiple of the oscillator frequency, it is nec- 
essary to take connections between different manifolds 
into account. Our approach is valid at resonance as well 
as positive and negative qubit detuning, provided that 
A < ri and/or strong coupling g/Q. 
As we already stated above, for weak coupling strengths 
like g/i} ~ 10~^ our approach cannot represent a replace- 
ment to the exactly solvable Jaynes-Cummings model. 
Also in the intermediate range of g/fl ^ 10~^, perturba- 
tivc approaches or the GRWA for zero bias calculation 
might be preferable. They fail, however, in the case of 
even stronger coupling - especially if the qubit is tuned 
away from its symmetry point. Here, our method shows 
that a new physical behavior can be expected, for which 
first hints have been given in recent experiments. 
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A^, 
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k^{j+l.j'+l} 
1 



-e + [k - j)n -e+{k- j')n 



Af+'Aj+' 



a^:+'a^:+' 



-e + {f + l-j)n ~e+{j+l-f)Q 



(A5) 



Using the above expressions, we find the eigenstates of 
H to second order in A as 



Ackno'wledginents 



i$±,,> = i<i>L°i> + i$S> + i$S> 



(A6) 
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Appendix A: Eigenstates obtained by Van Vleck 
perturbation theory 



with 



2 



|$W.)^sin^^|t,/>(t,J-'|i5(^'|;,J> 
.sign(Af ') cos^ ^ \i,f){i,j'\iS^''^\{J+l) 



j'=0 



(A7) 



For A = the eigenstates of the qubit-oscillator system 
read 



|t,J-)=ENgnO--/)]l^'-^'l4ML'}("/4)lt,/), 

(Al) 

iD) = ^[sign(/ - j)]|^'~^i4m],,'}(«/4)I ;,/>• 

(A2) 

The matrix elements of the Van Vleck transformation S 
in the basis of these states are to first order W6i : 



(j,;/t|i^(')|j',t/;>=- 



(±)b'-il A 



2 sTif-m 
x{l-6,±i,,,). (A3) 



and 

Ql °° _ 

|<^,)=sin^^i;,f)(;,j'|i5(2)|;,j-) 

+ Sign (Af ) cos ^ J2 It, j'Xt, j'li^'^'ltTTo 
j'=a 

10'°° - ~ 

-2'"'^ ^ i;,.7')a,/|i5(i)|t,fc')(t,A:'|i^(i)|;,j) 

j'=o,k'=a 

el- 



isign(Af' )cos^ 



j'=0,k'=0 



(A8) 



For |<&_|! ,) one just replaces sin -j- 



^+,3 



cos -J- and 
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Appendix B: Eigenenergies for e = using VVP 



We perform the summation in Eq. (J19p and show ana- 
lytical expressions for the first four energy levels obtained 
from VVP for the zero static bias case e = 0: 



E^.o =h 



■L 



A^c- 



40 



-(r(0,-a) + ln(-a)+7) 



T^|Ae-"/2| 



(Bl) 



In this representation, we neglect now the remote matrix 
elements, which turn out to yield fast rotating contribu- 
tions for A w 17 . A more elaborated justification is given 
in |27| . This procedure is quite similar to the standard 
rotating-wave approximation and we end up again with 
block-diagonal matrix. 



Eta 



n 



l + 7-Fe"(a-l) 



n 40 

a[a ~ 7(a - 2)] + (a - l)^[r(0, -a) + In(-a)] 



T-|A(l-a)e-"/2| 



(B2) 



where we used the Euler-Mascheroni constant 7 and the 
incomplete F- function plj . 



Appendix C: The generalized rotating-wave 
approximation (GRWA) 

Since we use the generalized rotating-wave approxima- 
tion in Sec. IIII Al where we calculate the energy spec- 
trum at £ = as a comparison to our VVP results, we 
will sketch its derivation in this appendix. A detailed 
description is found in [27[. The first step in its deriva- 
tion is to represent the qubit oscillator Hamiltonian ([1]) 
in the effective basis states (fT4|) and (fT5|) . disregarding 
the second-order corrections in A. Taking into account 

that A-? = (— l)l^~^ 'a-!,, the corresponding matrix is for 
the first six basis states {|$V ,)} with j = 0, 1, 2 
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(C2) 
which is straightforwardly diagonalized. The energy of 
the ground state remains unchanged, namely E-^o- The 
remaining levels are 



^GRWA/^ =(j + -)0 



r , A 



T-J 



e-"/2(m"(a)|-m"+i(")l) 



T 



1 -2 

)0-- , 
2^ 4 

|-^e-"/2(m"(«)i+mv(")i) 

A2 a 



-a ri 1 



4 j + l 



Lj(«)] 



(C3) 
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